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I . INTRODUCTION 

A. BACKGROUND 

Autonomous vehicles suitable for use in modern 
applications require high maneuverability, quick response, and 
robust performance characteristics [Ref. 7]. In order to 
maintain accurate track keeping characteristics, a suitable 
combination of path planning, navigation, guidance, and 
control design is needed [Ref. 10] . Sufficient information is 
obtained from charted obstacles and the environment for smooth 
paths to be generated for the vehicle to follow [Ref. 8]. 
Although it is possible to design a nonlinear controller which 
incorporates and utilizes information on the nonlinear dynamic 
properties of the vehicle, as well as the geometric 
nonlinearities of the desired track, the resulting scheme 
tends to be very complex and time consuming [Ref. 3]. The 
alternative is to use separated navigation, guidance, and 
autopilot functions. A certain level of feedback is provided 
at the navigation level through the use of sonars in order to 
replan a path due to uncharted obstacles or changes in mission 
objectives. This operation is event -driven and occurs 
asynchronously, and in many cases it does not affect stability 
and performance of the overall navigation, guidance, and 
control scheme. Based on the provided navigational 
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information, the guidance law generates heading commands, 
which are executed by the control law by appropriate use of 
vehicle effectors, such as control surfaces and cross body 
thrusters. However, the time required to process sonar data 
and/or inertial position information may be significant and 
cannot be neglected [Ref. 9] . In addition, the guidance and 
control laws must be as fast as possible in order to ensure 
accurate path keeping characteristics. The navigational 
position information time lags, as well as the dynamic lags 
due to the vehicle inertia, however, set a limit on the 
vehicle reaction time. Therefore, stability of the combined 
scheme becomes an issue which requires analysis. 

B. OBJECTIVE 

Previous studies have established boundaries for guidance 
and control laws in the horizontal [Ref. 10] and vertical 
plane [Ref. 12] maneuvering along straight line paths, as well 
as curved segments [Ref. 11]. The most important assumption 
in those results was the existence of instantaneous positional 
information updates when needed. In this work we relax the 
latter assumption. Stability boundaries are computed 
parametrized by the amount of positional information time lag. 
Results are presented based on eigenvalue and frequency 
response techniques. This thesis builds on the linear 
stability analysis performed in [Ref. 13] which can predict 
the shape of the stability boundaries. In this work, vehicle 
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motions in the vicinity of the corresponding boundary are 
assessed using nonlinear bifurcation theory [Ref. 5]. It is 
shown that the loss of stability occurs as generic Hopf 
bifurcations, where upon loss of stability of straight line 
motion a family of periodic solutions appears. Nonlinear 
expansions utilizing center manifold approximations and 
integral averaging techniques, reveal that these Hopf 
bifurcations can occur either in supercritical or subcritical 
forms. In the supercritical case, a stable family of limit 
cycles is generated immediately after the nominal motion 
becomes unstable. In the subcritical case, however, the 
resulting limit cycles are initially unstable and they are 
generated even before the nominal motion loses its stability. 
In the latter case, the domain of convergence of straight line 
motion becomes increasingly smaller as the stability boundary 
is reached. As a result, the methods developed in this work 
characterize the degree of confidence of the computed 
stability boundaries by examining stability under finite 
external disturbances. A first order memory scheme, which can 
be easily implemented in real time, is suggested in order to 
expand the region of stability of straight line motion. All 
computations are performed for the NFS autonomous underwater 
vehicle for which a set of geometric properties and slow 
motion hydrodynamic derivatives is available [Ref. 1] . Unless 
otherwise mentioned, all results are presented in standard 
nondimensional form. Equation development presented in 
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[Ref .13] has been used in chapter II, and with modifications 
in chapter IV. 

C. THESIS OUTLINE 

Chapter II presents the formulation of vehicle equations 
of motion and the criterion for determining the region of 
guidance/control stability in the presence of a position 
information time lag. 

Chapter III examines the vehicle's control stability 
through the use of a Hopf bifurcation analysis. 

Chapter IV examines the effect on the region of stability 
determination by the use of a memory model which incorporates 
the two previous vehicle positions. 
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II. PROBLEM FORMULATION 



A. EQUATIONS OF MOTION 

The mathematical model which describes vehicle motion in 
the horizontal plane consists of nonlinear sway and yaw 
equations of motion. A moving coordinate frame which is fixed 
at the vehicle's geometric center gives Newtonian equations 
of motion which are 

miv+ur+x^r-y^r^) =Y (2.1) 

iv+ur) -my^vi =N (2.2) 

where u is constant forward speed; v and r are the relative 
sway and yaw velocities of the moving vehicle relative to the 
water; m is the mass of the moving vehicle; Xq, y^ are the 
respective lateral and longitudinal positions of the center of 
gravity; and Y, N represent the total excitation sway force 
and yaw moment, respectively. The vehicle's added mass, 
damping, and viscous drag due to motion through the water are 
represented by quadratic drag terms and memoryless polynomials 
in V and r. Y and N can be represented as the sum of these 
terms by 
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Y= ^ ( Y^v+Y^ur) l^Y^uv 



N=-^l^Nfr+£-l^ (N^,v+N^ur) *£-l^N^uv 



where 1 is the vehicle length, Y^, represent partial 
derivatives of Y and N with respect to a, is the drag 
coefficient, and 6 is the rudder angle. 

At relatively high speeds and low angles of attack the 
steering response is predominantly linear. In low speed 
maneuvering or hovering conditions the cross flow integral 
drag term becomes significant. 

Vehicle yaw rate is expressed by 

ilr = r (2.3) 



and inertial position rate is expressed by 

y = usinilJ+vcosiir (2.4) 

where \J/ is the heading angle and shown in Figure 1. 
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Figure 1. Vehicle Geometry and Definition of Symbols. 
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Taking equations (2.1), (2.2), and (2.3) as a set of three 
coupled linear differential equations, and a linearized form 
of equation (2.4), the final set of equations of motion for 
steering control are 

= r (2 . 5a) 

v=aiiUv+ai2ur+jbiU^5 (2.5b) 

r = a 2 ^uv+a 22 ur+b 2 U-b (2.5c) 

y = uilr + v (2 . 5d) 



at any nominal forward speed. 

B. GUIDANCE AND CONTROL 
1. Guidance Scheme 

An orientation based command scheme is used where the 
commanded vehicle heading angle is a function of the line 
of sight angle a between the actual vehicle position and a 
reference position on the nominal path located at a constant 
lookahead distance d. Schematically, this is presented in 
Figure 1 . 

The line of sight angle o is defined by 
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tana = - — , 
d 



( 2 . 6 ) 



The autopilot will be called upon to orient the 
vehicle to the commanded heading angle which, in pursuit 

guidance, will equal the line of sight angle. This defines 
as 

il/^= -arctan^ {2,1) 



For relatively small angles 




( 2 . 8 ) 



2. Controller Design 

The steering control governing equations can be 
represented in the form 



x=Ax-^bb (2.9) 

where the state vector equation is 

x=[ij>,v,r]^ (2.10) 

Linear full state feedback is introduced in the form 
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6 = (ijr 



( 2 . 11 ) 



The gains k. 2 , and k 3 are computed by pole placement 



to give the desired dynamics to the closed loop system. The 
vehicle's longitudinal axis is pointed toward the desired 
heading by the difference in the control law. 

With a dimensionless time constant t^ the general form 
of the characteristic equation is 



X^+a3^X^ + a2^-^a3=0 



( 2 . 12 ) 



where 




1 



The controller gains are computed from 






(2.13a) 






(2.13b) 



b^u'^k^-'-b^u^k^ = -o-i~ ^ 



(2.13c) 
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These gain values are continuously updated due to changing 
nominal forward speeds. Figure 2 shows the three controller 
gains at unit forward speed versus the system time constant. 

C. TIME LAG 

All of the required state information for vehicle control 
is available at sufficient rates with the possible exception 
of the y position information. Due to time requirements for 
reduction of sonar returns and inertial navigation 
information, there may be a time lag in the y position 
information . 

This time lag T (sec) is introduced to the guidance 
control law and the commanded rudder angle becomes 



-arctan 



yit-T) 

d 



(2.14) 



For small angles 




(2.15) 



and the resulting linearized control law becomes 






(2.16) 
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CONTROLLER GAINS 




Figure 2 . Controller Gains versus at unit forward speed. 
(k]^=solid, k2=dashed, k3=dotted) 
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The linearized equations of motion become 



4r=r 



(2.17a) 



+ r+b^u^^yi t-T) 



(2.17b) 



i = b^u'^k{^-^ (a 2 iU+jb 2 U^ic 2 ) v 

+ (a^^u-^b^u^kj) r-^b^u^^yi t-T) 



y = ui[i + i^ 



(2.17d) 



D. STABILITY ANALYSIS 

Control law stability is determined by the eigenvalues of 
the matrix A from the linear system 

x = Ax (2.18) 



with the state vector 



x= [ijJ, v,r,y] ^ 



(2.19) 
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where the state equations have been linearized around straight 



line motion. The characteristic equation of (2.18) is a 
quartic of the form 



where the coefficients B, C, D, and E are functions of vehicle 
properties, controller gains, and lookahead distance. 

The characteristic equation will give a pair of complex 
conjugate roots which cross the imaginary axis under the 
following conditions. 



Minimum lookahead distance ^cnt stability can be 

determined by translating these conditions to 






( 2 . 20 ) 



BCD-B^E-D^ = 0 



(2.21a) 




(2.21b) 



=0 



(2.22a) 




(2.22b) 



where 
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(2.23a) 



a, =a^a,-a. 



(a^a2 {Jd^3.22 ^2^12 ^2^ 

“2 H il 









=*21 



-OL^U 



(2.23b) 



- (jbia22-‘^2‘3i2"^2) ^^1^22"^2^12~^2) ^3 

^•^2^11 ”■^1^21 ^ LZ 



(2.23c) 



Analysis of the system when operating with an information 
time lag requires approximation by either a first order Taylor 
expansion 

y{t-T) =y-Ty (2.24) 

a second order Taylor expansion 

y{t-T) =y-Ty+^y (2.25) 



a third order Taylor expansion 

rp2 rp3 

y(t-T) ^y-Ty*—y-—y 



(2.26) 



or a frequency response analysis using the Nyquist criterion. 
Figure 3 shows the region of stability given by each of these 
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methods with a one second (dimentionless) time lag. It can be 
seen that the agreement is, in general, very good among the 
three Taylor series approximations, especially as the time 
constant t^ is increased. The third order approximation 
coincides with the exact solution from the frequency response 
(Nyquist criterion) method. The first order method requires 
the highest value of d for stability at a given t^, and 
therefore is the most conservative method for design use. For 
this reason the first order method is used in the next chapter 
to study the dynamics of the system as the critical stability 
boundary is crossed. 
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d 

Figure 3. Taylor expansion and Nyquist stability analysis 
with a one second time lag. (Taylor first 
order=solid, Taylor second order=dashed, Taylor 
third order and Nyquist=dotted) 
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III. HOPF BIFURCATION 



A. INTRODUCTION 

It can be numerically confirmed that in all cases of 
stability loss of the previous chapter, one pair of complex 
conjugate eigenvalues of the corresponding eigenvalue problem 
crosses transversely the imaginary axis. A situation like 
this in which a certain parameter is varied such that the real 
part of one pair of complex conjugate eigenvalues of the 
linearized system matrix crosses zero, will result in the 
system leaving its steady state in an oscillatory manner. 
This loss of stability is called Hopf bifurcation and 
generically occurs in one of two ways, supercritical or 
subcritical. In the supercritical case, stable limit cycles 
are generated after the nominal straight line motion loses its 
stability. The amplitudes of these limit cycles are 
continuously increasing as the parameter distance from its 
critical value is increased. For small values of this 
criticality distance the resulting limit cycle is of small 
amplitude and differs little from the initial nominal state. 
In the subcritical case, however, unstable limit cycles are 
generated before the nominal state loses its stability. 
Therefore, depending on the initial conditions it is possible 
to diverge away from the nominal straight line path and 



18 



converge toward a different steady state even before the 
nominal motion loses its stability. In many cases this new 
steady state is another limit cycle of considerably higher 
amplitude. This means that in the subcritical Hopf 
bifurcation case the domain of attraction of the nominal state 
is decreasing and in fact it shrinks to zero as the critical 
point is approached. Random external disturbances of 
sufficient magnitude can throw the vehicle off to an 
oscillatory steady state even though the nominal state may 
still remain stable. After the nominal state becomes 
unstable, a discontinuous increase in the magnitude of motions 
is observed as there exists no simple stable nearby attractors 
for the vehicle to converge to. Distinction between these two 
qualitatively different types of bifurcation is, therefore, 
essential in design. The computational procedure requires 
higher order approximations in the equations of motion and is 
the subject of the next section. 

B. COMPUTATIONS 

The vehicle equations of motion are written in the form 

i|f =r (3 . la) 

v = aj^^uv-^a^2^r-^b^u^b (3.1b) 
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(3.1C) 



y = usinilF+vcosi|r 



(3. Id) 



In order to capture the effect of rudder saturation we 
write the rudder angle 6 in the form 



where is the saturation limit of 6 , typically set at +0.4 
radians. Equation (3.2) is used instead of a hard saturation 
function because of its smoothness, which is important for the 
computation in this section. 6 q is the linear control law of 
the form 



using the first order approximation for y(t-T). 

The state equations (3.1) through (3.4) are written in the 
compact form 



6=6^^, tanh(-^) 

^ sat 



(3.2) 






(3.3) 



where 




(3.4) 
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x=f{x) 



(3.5) 



where f (x) is a nonlinear function of the state variables 
vector 

, V, I , y] (3.6) 

Expanding f (x) in Taylor series, we can write (3.5) in the 
form 



x=Ax+p(x) (3.7) 

where A is the linear system matrix and g (x) contains all the 
leading nonlinear terms of f (x) . Due to the port/starboard 
symmetry of our problem, g (x) contains only third order terms. 
Defining T as the matrix of eigenvalues of A and applying the 
transformation 

x=Tz (3.8) 



the state equation is transformed into its canonical form 

z = T-^ATz + T-^g (Tz) (3.9) 



At the Hopf bifurcation point 



21 



0 -(i)« 0 0 




(3.10) 



0 0 p 0 

0 0 0 g 



where p and q are negative and coq is positive. In the new 
system of coordinates 



the dynamics of the system are generated by a reduced two 
dimensional system 7 .^ and 7 . 2 , since the coordinates Z 3 and 
correspond to the eigenvalues p and q and are asymptotically 
stable. These stable variables Z 3 , z^ can be expressed as a 
function of the critical variables 7 -^, Z 2 ,* these functional 

expressions are at least of third order. Therefore, the 

stable coordinates Z 3 , z^ do not influence the third order 

expansion of (3.7). Using the above expression, we can write 

the two dimensional system in z^, Z 2 as 



where the coefficients r^^j are complicated expressions that 
are derived from ( 3 . 9 ) 



(3.11) 




(3.12a) 




(3.12b) 
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Equations (3.12) are valid at exactly the Hopf 
bifurcation point. For values of the parameter d close to the 



bifurcation point, they are written as 



- (G)Q-^-G)'e) {z^, z^) 



(3.13a) 



^2 ” Z^-^Ol'zZ2-^F2 {z^, Z^) 



(3.13b) 



where of oj ' are the derivatives of the real and imaginary 
parts of the critical pair of eigenvalues evaluated at the 
bifurcation point; ^ is the difference of the parameter d from 
its critical value; and the nonlinear function , F 2 remain 
the same as in (3.12), 





If we introduce polar coordinates in the form 



=i?COS0 



(3.15a) 



Z 2 =i^sin0 



(3.15b) 



equations (3.13) become 
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R = a^tR+F^ {R, 0) COS 0 +F 2 (R, 0) sin0 



(3.16a) 



= ((Oq+w'e) F+F 2 (R, 0) cos0-F^ {R, 0) sin0 (3 . 16b) 



Equation (3.16a) yields 



R = a'tR+P(d) R^ 



(3.17) 



where P ( 0+27 t) =P ( 0 ) . If (3.17) is averaged over one cycle in 
6, we get an equation with constant coefficients 



Carrying out the indicated integration in (3.19) we obtain 



Existence and stability of limit cycles can be determined 
by analyzing the equilibrium points of the averaged equation 
(3.18), which correspond to periodic solutions in z^, Z 2 as 
can be seen from (3.15) . If we examine equation (3.18) we can 
see that: 



R = a'eR+KR^ 



(3.18) 



where 




(3.19) 




(3.20) 
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1) If O' '>0, then: 

(a) If K>0 then unstable limit cycles coexist with the 
stable equilibrium for £<0. 

(b) If K<0 then stable limit cycles coexist with the 
unstable equilibrium for £>0. 

2) If Ckf '<0 , then: 

(a) If K>0 then unstable limit cycles coexist with the 
stable equilibrium for £>0. 

(b) If K<0 then stable limit cycles coexist with the 
unstable equilibrium for £<0. 

Therefore, we can see that computation of the above 
nonlinear coefficient K can distinguish between the two 
different types of Hopf bifurcation: 

• Supercritical if K<0 

• Subcritical if K>0 
[Ref. 10]. 

C . RESULTS 

Values of the coefficient K have been calculated for 
several different values of a y position information time lag, 
over a span of system time constants, using the Fortran 
program presented in appendix A. This program has also been 
used to determine the vehicle period of oscillation in the 
same conditions. Figure 4 shows the behavior of the K 
coefficient at three different time lags, and Figure 5 shows 
periods of oscillation at the same three time lags. It can be 
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Figure 4. Coefficient K versus at three values of time lag 
T. (0.5 secasolid, 1.0 sec=dashed, 1.5 sec=dotted) 
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PERIOD 




tc 



Figure 5. Vehicle oscillation period versus for three 

values of time lag T. (0.5 sec=solid, 1.0 
sec=dashed, 1.5 sec=dotted) 
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seen that supercritical bifurcations are ensured for 
sufficiently high values of the control time constant t^. For 
t^ less than a certain critical value, the corresponding Hopf 
bifurcations are subcritical. This situation should be 
avoided in practice. 

D. SIMULATIONS 

The vehicle path has been simulated by using the Euler 
integration method coded in a Fortran program presented in 
appendix B. The vehicle control law (2.17) as well as the 
control gains (2.13) are used. The program has been run with 
input values for system time constant, y position information 
time lag, and lookahead distance. The vehicle was given an 
initial nominal y offset of half a shiplength and a nominal 
forward speed. The resulting plots of y position against time 
show the vehicle's oscillatory path either converge to the 
commanded path or diverge. 

The results of these simulation runs were compared to the 
results of the Hopf bifurcation computations in two ways. The 
period of oscillation observed in the vehicle path was 
compared to the period predicted by the Hopf program. 
Additionally, the vehicle stability, determined by convergence 
to the commanded path, was compared to the K coefficient sign 
prediction of stability which was determined by the Hopf 
computations. Figure 6 shows a vehicle path in stable 
conditions. Figure 7 shows a vehicle path in unstable 
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conditions . 



These simulations were chosen 



in the 



supercritical bifurcation case and it can be seen that 
period of oscillation observed from Figures 7 agrees very 
with the theoretical results of Figure 5. 



the 

well 
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I 

Figure 6. Vehicle path in stable conditions where time 
constant t^, = 1.0 sec, time lag T = 1.0 sec, and 
lookahead distance d = 2 shiplengths. 
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t 

Figure 7. Vehicle path in unstable conditions where time 
constant t^ = 1.0 sec, time lag T = 1.0 sec, and 
lookahead distance d = 1.25 shiplength. 
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IV. 



STABILITY ENHANCEMENT 



A. TWO -POINT FORMULA 

The time lag of y position information discussed in 
chapter II has thus far been represented by a single piece of 
inf onnation used in the control law with some delay assigned 
to it. In an attempt to improve stability with regard to this 
delay, the use of the previous two positions in the control 
was examined. 

Equation (2.15) shows the representation of the time 
delay, and equation (2.16) shows its effect on the control 
law. A straight line approximation is applied to the previous 
two positions as shown in Figure 8. The general y position is 
defined by 



y ( t) =a^t+a^ 



(4.1) 



The two previous positions become 






(4.2a) 






(4.2b) 



The coefficients and aj are stated as 
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Figure 8. 



Straight line approximation using 



two previous 



positions . 






^ 1-^2 



(4.3a) 



. _ yi^2~y2^i 

^2 TTT 



(4.3b) 



The previous time values are represented as 



(4.4a) 



t^ = t-T (4.4b) 

where T is the time lag associated with the y position 
information. The previous positions, in terms of the previous 
time values are 



y^=y{t-2T) 



(4.5a) 



y2=y( t-D 



(4.5b) 



Substitution and algebra gives a general position 
expression of 



y( t) =2y{t-T) -y{t-2T) 



(4.6) 



34 



B. STABILITY ANALYSIS 



The Laplace transform of equation (4.6) is 

y(t) =y(2e-’’^-e-2^®) <4.7) 

The control law becomes 

6 =JCi<|;+Jt2V+J^3r+-^y(2e-’’®-e-2’’^) (4.8) 

and the linearized equations of motion become 

i|j=r (4.9a) 

k (4.9b) 

+bj^u^-^y (2e'^^-e'^^®) 

r = b 2 U^k^}\!^ {a^ 2 ^-^b 2 U^k^) v-^ {a 22 U+b 2 U^k^) r 

^ > Tc 9TCV (4.9c) 

^ d 

y=L 2 t|j + v (4.9d) 

These motion equations reduce to the state space form 

x=Ax (4.10) 
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The matrix form becomes 



'f 




0 


0 1 0 




V 




b^u^k^ 


aiiU+jbiU^ic2 a^2^-^b^u^k2 


V 


f 




b^u^k^ 


^3 U2 ^ ( 2 e e 


r 


y 




u 


10 0 


y 



The characteristic equation of the form [A-Is]=0 is 

D{2e~'^^-e~^'^^) (4.11) 



where 

A 3 = - (< 311 +^ 22 ) u - ib^k^+b^k^) 



^2 = ^3i1^22-^12^2i'I <'3iA“'32A' “ ^■'^3 < ' 322*1 "' 312*2 > 



A = -( 321 *i"' 3 i 1*2) U '*1 



B2 = -jb^u^iCi 



Bi = -*2U^;Ci-('3i2*2 -332*1) “^*1 
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■®0 ^^ 11-^2 ^ 21 -^ 1 ^ 




The characteristic equation of the system is written as 

1+KGis) =0 (4.12) 



where 



s is^+A^s'^+A^s+A^) 



(4.13) 



is the transfer function, and we have denoted 



K=D 



(4.14) 



With s = jcj, the phase angle is represented by 

4) = ZG(ja)) (4.15) 



where 



-Z ( -ja>^-^3(i)^+A2ja)+^i) 
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^ 2cos(a)D -cos(2a)D 2 S.-B^co^ 



-tan'^ ( - 



-^) 



(4.16) 



-A^(i>^+A^ 



The Nyquist stability criterion states that at the solution of 

4> (o)) = -71 (4.17) 

where cj is the phase cross over frequency the gain margin 
must equal 1 



\KGij(D^) I =1 



(4.18) 



where K = D = 1/d. Equation (4.17) becomes 
, -2sin (o). D +sin (2g). T) 

T I ± ± 



tan -^tan 

2cos (o)^r) -cos (2o)^D 2 






-Ul+AjOJi 

-can r -= — -) =-n 

-Ajto^+Aj 



(4.19) 



Rearranged as 



tan ^ ( 



B,-B^<^1 



-tan"^ ( 



TC 

2 



•* 2 , 

-Aj (0 1 +A^ 

-2sin(co,r) +sin{2(j>.T) 

-tan ^ ^ : i— 

2cos (co^D -cos (2o)^D 
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Taking the tangent of both sides and using the identity 



tan {x-y) 



_ tan (x) -tan (y) 
1+tan ix) tan (y) 



then algebra and rearranging gives 

( 4 . 20 ) 

2cos (o)^r-cos (2(*)^D 
-2sin(a)iD +sin (2(*)j^D 

The stability computations require an initial value for 
This is accomplished by setting the time lag T=0, and 

using 

{B^-B^iSil) (-^30)1+^^) (-a)i+A2a)i) =0 ( 4 . 21 ) 

rearranged as 

(B2A3-B1) (B3A2-B2A3-B0A3) 0)?+B„Ai = 0 ( 4 . 22 ) 

Setting the gain margin equal to 1 gives 
1 |(-B3 0)^-^B3J0)3^Bq) (2e'^'‘-’'^-e~^^“-^) | 

and d can be solved for by 
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(4.23) 



^ (J9q-jB2 0 )?) +B^a)i y^5-4cos (co^T 

(i4^“i43G)i) (i42C03^-0)i) ^ 

In order to facilitate computation, terms are grouped as 
follows 



a^= B^A^ -B^A^ ~BqA^ 



a,=BaA, 



Pl = -^2 

P 2 ='®0^^2^2~'®A 



P,=b,a,-b,a, 



Equation (4.20) becomes 

Pia)^ + P 2 Ct)^ + P 3 a) _ 2cos (cjT) -cos (2(oT) 
+a 2 (t)^+a^ -2sin (wT) +sin (2a)!T) 



Put in the form 
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f (g>) =0 



(4.25) 



equation (4.24) becomes 

[ ( Pj^G)^ + P 20 )^ + P 3 (i)) (-2sin((i)D +sin(2(i)r) ) ] 

- [ (a^G)‘*+a 20 )^+a 3 ) (2cos(o>T) -cos (2o>T) ) ] =0 



Newton's iteration is used as 









Where the function of cj is 

f((D) = ( P^g>^ + P2(i>^ + P3W) (-2sin(o>r) ) 
+ (P^Ci)^ + P2(i)^-*-P30>) (sin(2a)D) 
- (a^o)^+a20)^+cc2^ (2cos(o>r)) 

+ (a^oy^+a^Ciy^-^a^) (cos(2a>D) 



and the derivative function is 

f^{(jy) = (5p^(i)^+3p2W^ + P 3 ) (-2sin(co7’) ) 

- (p^o)^ + P2a>^-^P^a)) {2TCOS {(dT) 

+ (5P^g)^+ 3P20>^ + P^) (sin (2 (i)7’) ) 

-t- ( P^G>^ + P20)^ + p3<*)) {2TCOS {2(jyT) ) 
+ (4a^(i)^+2a2(i>) ( -2cos (wT) ) 

+ (a^o)^+a20)^ + a3) { 2 Tsin{(jyT) ) 

+ (4ai(i)^+2a2(i)) (cos(2a)T)) 

- ia^(jy^+a2(jy^^a^) { 2 Tsin{ 2 uyT) ) 



(4.26) 



(4.27) 



(4.28) 



(4.29) 
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Appendix C presents the Fortran program used to determine 
stability characteristics as a function of system time 
constant and time lag. 

C. RESULTS AND DISCUSSION 

Values for the minimum lookahead distance required for 
controller stability have been computed as a function of the 
system time constant and the y position information time lag 
using the two-point formula. Figures 9, 10, and 11 show 
critical lookahead distance plotted against time constant at 
three values of time lag for both the two-point method as well 
as the single point (Nyquist) method discussed in chapter II. 
These figures show an overall decrease in required lookahead 
distance as a result of the use of the two-point memory model. 
An exception to this occurs at low values of t^. However, 
these values are to be avoided in practice, as explained in 
the previous chapter. 

Figure 12 through 15 show the critical lookahead distance 
plotted against time lag for different values of the time 
constant using both the two-point and single point methods. 
Again the two-point method produces superior results to the 
single point method, with the exception of small time 
constants and large time lags. The same information is 
summarized in Figure 16 through 19, where we show the ratio 
d 2 /di (critical lookahead distance using the two-point method 
divided by the critical lookahead distance using the single 
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point method) verses time constant and time lag. It can be 



seen that 
improvement 



the use of the two-point method can result in 
in the region of stability of over 20%. 
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d 

Figure 9. Critical lookahead distance d versus time constant 
t^ with time lag T = 0.5 sec (two-point 

methodssolid, single point=dashed) . 
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d 

Figure 10. Critical lookahead distance d versus time constant 
t^, with time lag T = 1.0 sec (two-point 

method=solid, single point=dashed) . 
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Figure 11. Critical lookahead distance d versus time constant 
with time lag T = 1.5 sec (two -point 

method=solid, single point^dashed) . 
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d 

Figure 12. Critical lookahead distance d versus time lag T 
with time constant t^ = 0.8 sec (two-point 
methodssolid, single points dashed) . 
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d 

Figure 13. Critical lookahead distance d versus time lag T 
with time constant t^ = 1.0 sec (two-point 

method=solid, single points dashed) . 
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Figure 14. Critical lookahead distance d versus time lag T 
with time constant t^ = 1.5 sec (two -point 

method=solid, single point=dashed) . 
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Figure 15. Critical lookahead distance d versus time lag T 
with time constant t^ = 2.5 sec (two-point 

method=solid, single point=dashed) . 
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Figure 16. Ratio ^ 2/^1 versus time constant t^^ with time 
lag = 1.0 sec. 
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Figure 17. Ratio versus time constant t^ with time 

lag = 1.5 sec. 
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figure 18. Ratio 6 . 2 / versus time lag T with time constant 
1.0 sec 
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Figure 19. Ratio versus time lag T with time 

constant = 1.5 sec. 
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V. 



CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

Nonlinear bifurcation theory has been applied to the NPS 
autonomous underwater vehicle through a Hopf bifurcation 
analysis of the vehicle's navigation guidance and control 
scheme. It has been shown that under normal conditions the 
vehicle will operate in the supercritical region, and 
therefore the unpredictable behavior associated with the 
subcritical case should not to be expected. 

Secondly, the control scheme was modified by the use of a 
single stage memory model which incorporates the two previous 
vehicle positions in the linearized control law. It has been 
shown that the use of this memory model reduces the vehicle's 
required lookhead distance for control stability within the 
normal operating range. 

B. RECOMMENDATIONS 

In the event that not all of the state information 
required is directly available through measurements, the 
control/guidance scheme would include an observer. The effect 
of the observer dynamics on the higher order nonlinear terms, 
and therefore on the bifurcation limit cycle stability, should 
be examined. 
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APPENDIX A 



PROGRAM HOPF.FOR 

Hopf Bifurcation Analysis 

Critical Point: Exact 

Third Order Expansions: First Order Approximation 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION Kl , K2 , K3 , K4 , L , NR , NV , NDRS , NDRB , I Z , MASS , 

& NRDOT,NVDOT, KlP, K2P 

DOUBLE PRECISION Ml 1 , Ml 2 , Ml 3 , Ml 4 , M2 1 , M22 , M2 3 , M2 4 , 

1 M31 ,M32 ,M33 ,M34 ,M41 ,M42 ,M43 ,M44 , 

2 N11,N12,N13,N14,N21,N22,N23,N24, 

3 N31,N32,N33,N34,N41,N42,N43,N44, 

4 L21,L22,L23,L24,L31,L32,L33,L34, 

5 L41,L42,L43,L44 

DIMENSION A(4,4),T(4,4) ,TINV( 4,4), FVl ( 4 ) , IVl (4),ZZZ(4,4) 
DIMENSION WR(4) ,WI(4) ,TSAVE(4,4) ,TLUD(4,4) ,IVLUD(4) ,SVLUD(4) 
DIMENSION ASAVE(4,4) 

OPEN ( 11 , FILE=' HOPF. DAT' , STATUS= ' NEW ' ) 

OPEN (12,FILE='PERI.DAT' , STATUS= ' NEW ' ) 

OPEN ( 13 , FILE='ALPH10 .DAT' , STATUS= ' NEW ' ) 

Vehicle Parameters 

WEIGHT=435.0 
IZ =45.0 

L =7.3 

RHO =1.94 

G =32.2 

XG =0.0104 
MASS =WEIGHT/G 
U =5.0 

RATIO =0.0 

DO =0.4 

WRITE (*,1006) 

READ (*,*) DO 

YRDOT =-0. 00178*0. 5*RHO*L**4 
YVDOT =-0. 03430*0. 5*RHO*L**3 
YR =+0 .01187*0 . 5*RHO*L** 3 
YV =-0 . 03896*0 . 5*RHO*L** 2 
YDRS =+0.02345*0. 5*RHO*L**2 
YDRB =+0.02345*0. 5*RHO*L**2 
NRDOT =-0 . 00047*0 . 5*RHO*L**5 
NVDOT =-0 . 00178*0 . 5*RHO*L**4 
NR =-0 . 01022*0 . 5*RHO*L**4 
NV =-0 . 00769*0 . 5*RHO*L**3 
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!j)RS J 

^)RB -+0.28 3*0.02 3 4 5*0.5*RHO*L**3 

:i -( IZ-NRDOT) *(WASS-YVDOT)- 

(MASS*XG-YRDOT) * (MASS*XG-NVDOT) 
till=( ( IZ-NRDOT) *YV-(MASS*XG-YRDOT) *NV)/DH 
ta2=( ( IZ-NRDOT) * (-WASS+YR)- 

( MASS*XG-YRDOT ) * ( -MASS*XG+NR ) )/DH 
Ia21=( (MASS-YVDOT) *NV-(MASS*XG-NVD0T) *YV)/DH 
U22=*( (MASS-YVDOT) * (-MASS*XG+NR)- 

( MASS*XG-NVDOT ) * ( -MASS+YR ) ) /DH 
|ill=( ( IZ-NRDOT) *YDRS-(MASS*XG-YRDOT) *NDRS)/DH 
[112* ( ( I Z-NRDOT ) * YDRB- ( MASS*XG-YRDOT ) *NDRB )/DH 
[121* ( ( MASS-YVDOT ) *NDRS- ( MASS*XG-NVDOT ) * YDRS )/DH 
[122- ( (MASS-YVDOT) *NDRB-( MASS *XG-NVDOT) * YDRB) /DH 



[11-BB11 + RATI0*BB12 
[12 = BB21 + RATI0*BB22 

•UTE (*,1005) 

[SAD ( * , * ) TL 
'UTE (*,1001) 

•SAD ( * , * ) TC 



>S*1 .D-10 
[?MAX=10 00 

PL*0.5 

>0 50 M=l,100 
■.D=TL*L/U 

:c=o.5 

) 40 K*l,100 
TCD*TC*L/U 
POLEC-l.O/TCD 
AD1=3.0*POLEC 
AD2-3.0*POLEC**2 
AD3*POLEC**3 
A1*BB1*U*U 
B1=BB2*U*U 

Cl=-ADl-(AAll+AA22 ) *U 
A2=(BB1*AA22-BB2*AA12 ) *U**3 
B2=(BB2*AA11-BB1*AA21 ) *U**3 
Kl*AD3/( ( BB2*AA11-BB1*AA21 ) *U**3 ) 
C2*AD2-(AA11*AA22-AA12*AA21 ) *U* * 2+BB2 *U*U*Kl 
K2*(C1*B2-C2*B1 )/(Al*B2-A2*Bl ) 

K3=(C2*A1-C1*A2 )/( A1*B2-A2*B1 ) 

Al* (AA11*BB2-AA21*BB1)*U*U*U*K1 

A2* (AA11*AA22-AA12*AA21 ) *U*U+ ( AAl 1 *BB2-AA2 1 * BBl )*U*U*U*K3 + 
(AA22*BB1-AA12*BB2 ) *U*U*U*K2-BB2*U*U*K1 
A3*-(AAll+AA22 ) *U- ( BBl *K2+BB2*K3 ) *U*U 
BO* (AA11*BB2-AA21*BB1 ) *U*U*U*U*K1 
B1*-(AA12*BB2-7^22*BB1 ) *U*U*U*K1-BB2*U*U*U*K1 
B2—BB1*U*U*K1 

Compute Initial Approximation For Omega (TL=0) 
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AQ=b2*A3-B1 

BQ«B1*A2-B2*A1-B0*A3 

CQ=BO*Al 

DET=BQ**2-4 .D0*AQ*CQ 
IF (DET.LT.O.DO) GO TO 40 
ZTl«(-BQ+DSQRT(DET) )/(2.0D0*AQ) 

ZT2«( -BQ-DSQRT( DET) )/( 2 . 0D0*AQ) 

IF ( ZTl .GE. 0 .DO ) OM=DSQRT( ZTl ) 

IF ( ZT2 .GE . 0 .DO ) OM«DSQRT( ZT2 ) 

ALPHAl-AQ 
( ALPHA2=BQ 

( ALPHA3«CQ 

{ BETAl«-B2 

( BETA2-B0+B2*A2-B1*A3 

( BETA3*B1*A1-B0*A2 

r 

r Compute Exact Omega Using Newton's Iteration 



93 

& 

& 

& 

& 

& 



C 
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F=( BETA1*0M0LD**5+BETA2*0M0LD**3+BETA3*0M0LD) *dsin( omold*tld ) 
+ ( ALPHAl *OMOLD* * 4+ALPHA2 *OMOLD* * 2+ALPHA3 ) *DCOS ( OMOLD*TLD ) 
FPRIME=( 5 .DO*BETA1*OMOLD**4+3 . DO*BETA2*OMOLD* * 2+BETA3 ) 
*DSIN(OMOLD*TLD)+ ( BETAl *OMOLD* * 5+BETA2 *OMOLD* * 3+BETA3 *OMOLD ) 

*TLD*DCOS(OMOLD*TLD)+( 4 .DO*ALPHAl*OMOLD**3+2 .DO*ALPHA2*OMOLD 

*DCOS ( OMOLD*TLD ) - ( ALPHAl *OMOLD* * 4+ALPHA2*OMOLD* * 2+ALPHA3 ) 
*TLD*DSIN( OMOLD*TLD) 

IF ( FPRIME . EQ . 0 . DO ) STOP 1112 
IF (F.EQ.O.DO) GO TO 92 
OMNEW=OMOLD-F/FPRIME 
OMDI F=DABS ( OMNEW-OMOLD ) 

IF (OMDIF.LE.EPS) GO TO 92 

OMOLD=OMNEW 

IT=IT+1 

IF ( IT.GT. ITMAX) STOP 1111 
GO TO 93 



) 



OM»OMNEW 

XDNUM*DSQRT( ( B0-B2 *OM* * 2 ) * * 2+ ( Bl *OM ) * * 2 ) 
XDDEN=OM*DSQRT( ( A1-A3 *OM* * 2 ) * * 2+ ( A2 *OM-ON* * 3 ) **2 ) 
XD=XDNUM/XDDEN 



Start Hopf Bifurcation Analysis 

Evaluate Nonlinear Rudder Expansion Coefficients 



C 



KlP=Kl-Kl*TLD*U/XD 

K2P=K2-K1*TLD/XD 



& 

& 



& 

& 



DPPV*- ( 1 . D0/( 3 . D0*D0 * * 2 ) ) * 3 . DO* K1P*K1P* K2P 

+ 0 5D0*Kl*TLD/XD + 3 . DO* Kl * ( TLD*U ) * * 3/ ( 3 . D0*XD* * 3 ) 
DPVV=- (i.D0/( 3.D0*D0**2) ) *3 .D0*K1P*K2 P*k 2P 
+ 3 . DO*U*TLD*TLD*TLD*Kl/( 3 . D0*XD** 3 ) 

DPPR=- ( 1 . D0/( 3 .D0*D0**2 ) ) * 3 . D0*KlP*KlP*K3 
DPRR=- ( 1 .D0/( 3 .D0*D0**2) ) *3 .D0*K1P*K3*K3 
DPPY=- (1.D0/(3.D0*D0**2) )*3. D0*KlP*KlP*Kl/XD 
- 3 D0*TLD*TLD*U*U*Kl/( 3 .D0*XD**3 ) 

DPYY*- (l.D0/( 3.D0*D0**2) ) * 3 . D0*KlP*Kl *Kl/ ( XD* * 2 ) 

+ 3 .D0*TLD*U*Kl/( 3 .D0*XD**3 ) 

DVVR=- ( 1 .D0/( 3.D0*D0**2) )*3.D0*K2P*K2P*K3 
DVRR=- ( 1 .D0/( 3.D0*D0**2))*3.D0*K2P*K3*K3 
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- 3.D0*K1*TLD*TLD/(3.D0*XD**3) 

DVYY=- ( 1 .D0/( 3 .D0*D0**2 ) ) * 3 . DO *K1 * Kl *K2P/ ( XD* * 2 ) 

+ 3.DO*TLD*Kl/( 3.D0*XD**3) 

DRRY=- ( 1 .D0/( 3 .D0*D0**2 ) ) * 3 . DO *K1 *K3*K3/XD 
DRYY=- (l.DO/( 3.D0*D0**2) ) * 3 . DO *K1 * Kl * K3/( XD* * 2 ) 

DPVR=“ ( 1 .DO/( 3 .D0*D0**2 ) ) *6 .D0*K1P*K2P*K3 
DPVY»- ( 1 .DO/( 3 .D0*D0**2 ) ) *6 . DO*KlP*Kl *K2P/XD 

- 6 .DO*TLD*TLD*U*Kl/( 3 .DO*XD**3 ) 

DPRY=- ( 1 .DO/( 3 .D0*D0**2 ) ) *6 . DO*KlP*Kl *K3/XD 
DVRY*- ( 1 .DO/( 3 ,D0*D0**2 ) ) *6 . DO*Kl *K2P*K3/XD 
DPPP=- ( 1 .DO/( 3 ,D0*D0**2 ) ) *1 . DO*KlP*KlP*KlP 

+ Kl*TLD*U/(6.D0*XD) + ( Kl * ( TLD*U ) * * 3 )/( 3 .D0*XD**3) 
DVW=- ( 1 .DO/( 3 .D0*D0**2 ) ) *1 . D0*K2P*K2P*K2P 
+ K1*(TLD**3)/(3.D0*XD**3) 

DRRR=- (l.DO/(3.DO*DO**2) ) *1 . D0*K3*K3*K3 

DYYY=- (1 ,DO/( 3.D0*D0**2) ) * ( Kl/XD ) * * 3-Kl/( 3 . DO *XD* * 3 ) 

- Kl/(3.D0*XD**3) 

Evaluate Transformation Matrix of Eigenvectors 

PSI =0.D0 

Y =0.D0 

V =0.D0 
DPSI=KlP 
DV =K2P 
DR =K3 

DY =K1*XD/(XD**2+Y**2 ) 

A(1,1)=0.0D0 
A(1,2)=0.0D0 
A(1,3)=1.0D0 
A( 1 , 4 )=0 . ODO 

A(2,1)= BB1*U*U*DPSI 

A(2,2)=AA11*U+BB1*U*U*DV 
A(2,3)=AA12*U+BB1*U*U*DR 
A(2,4)= BB1*U*U*DY 

A(3,1)= BB2*U*U*DPSI 

A( 3, 2 )=AA21*U+BB2*U*U*DV 
A( 3, 3 )=AA22*U+BB2*U*U*DR 
A(3,4)= BB2*U*U*DY 

A(4,l )=U*DCOS(PSI )-V*DSIN(PSI ) 

A(4,2)= DCOS(PSI) 

A(4,3)=0.0D0 
A(4,4)=0.0D0 
DO 11 1=1,4 
DO 12 J=l,4 

ASAVE( I , J)=A( I , J) 

CONTINUE 

CONTINUE 

CALL RG(4,4,A,WR,WI,1,ZZZ,IV1, FVl , I ERR ) 

CALL DSOMEG( IEV,WR,WI , OMEGA, CHECK) 

OMEGAO=OMEGA 
DO 5 1=1,4 

T(I,1)= ZZZ(I,IEV) 

T( I , 2 )=-ZZZ( I , IEV+1 ) 

CONTINUE 

IF ( lEV.EQ. 1 ) GO TO 13 

IF (IEV.EQ.2) GO TO 14 

IF (IEV.EQ.3) GO TO 15 

STOP 3004 
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T(I,3)=ZZZ(I,1) 

T( 1 , 4 )=ZZZ( 1 , 4 ) 

6 CONTINUE 
GO TO 17 

15 DO 7 1 = 1,4 

T(I,3)=ZZZ(I,1) 

T(I,4)-ZZZ(I,2) 

7 CONTINUE 
GO TO 17 

13 DO 16 1=1,4 

T(I,3)=ZZZ(I,3) 

T(I,4)=ZZZ(I,4) 

16 CONTINUE 

17 CONTINUE 

Normalization of the Critical Eigenvector 
CALL NORMAL(T) 



( 

{ 

< 

( 



C 



Definition of Mij 

M11=T(1,1) 
M21=T(2,1 ) 

N31=T( 3,1 ) 
M41=T(4,1) 
M12=T(1,2) 
M22=T(2,2) 

M32=T( 3,2) 

M42=T( 4,2) 
M13=T(1,3) 
M23=T(2,3) 

M33=T( 3,3) 

N43=T( 4 , 3 ) 
M14=T(1,4) 
N24=T(2,4) 
M34=T(3,4) 

M44=T( 4,4) 

Definition of Lij 



L21= DPPV*Mll*Mll*M21 + DPVV*N11*M21*M21 + 
& + DPRR*Mll*M31*M31 + DPPY*Mll *Nll*M41 + 
& + DVVR*N21*M21*M31 + DVRR*M21 *M31 *M31 + 
& + DVYY*M21*M41*M41 + DRRY*N31 *M31 *M4 1 + 
& + DPVR*Mll*M21*M31 + DPVY*Ml 1 *M2 1 *M4 1 + 
& + DVRY*M21*M41*M31 + DPPP*Ml 1 *Ml 1 *Ml 1 + 
& + DRRR*M31*M31*M31 + DYYY*M4 1 *M4 1 *M4 1 



PPV = Mll*Mll*M22 
PVV = N12*M21*M21 
PPR = Mll*Mll*M32 
PRR = M12*M31*M31 
PPY = Nll*Mll*M42 
PYY = M41*M41*M12 
VVR = M21*M21*M32 
VRR = M22*M31*M31 
VVY = M21*M21*M42 
VYY = M22*M41*M41 
RRY = M31*M31*N42 



+ 2.0*M11*N12*N21 
+ 2. 0*Mll*M21*M22 
+ 2 . 0*N11*M12*M31 
+ 2 . 0*Mll*M31*M32 
+ 2.0*M11*M12*M41 
+ 2.0*M11*M41*M42 
+ 2.0*M31*M21*M22 
+ 2 . 0*M31*M32*M21 
+ 2 . 0*M41*M21*M22 
+ 2 . 0*M41*M42*M21 
+ 2 . 0*M41*M31*M32 



DPPR*Mll*Mll*M31 

DPYY*Mll*M41*M41 

DVVY*M21*M21*M41 

DRYY*M31*M41*N41 

DPRY*Mll*M41*M31 

DVVV*M21*M21*M21 
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RYY « + 2 . 0 *M4 1 *M4 2 *M3 1 

PVR = M11*M21*M32+M31*(M11*M22+M12*M21 ) 
PVY = M11*M21*M42+M41*(M11*M22+M12*M21 ) 
PRY = Mll*M41*M32+M31*(Mll*M42+Ml2*M41 ) 
VRY = M21*M41*M32+M31*(M21*M42+M22*M41 ) 
PPP - 3 . 0*M11*M11*M12 
VVV - 3.0*M21*M21*M22 
RRR - 3 . 0*M31*M31*M32 
YYY - 3 . 0*M41*M41*M42 



L22=DPPV*PPV+DPVV*PVV+DPPR*PPR+DPRR*PRR+DPPY*PPY+DPYY*PYY 

+DVVR*VVR+DVRR*VRR+DVVY*VVY+DVYY*VYY+DRRY*RRY+DRYY*RYY 

+DPVR*PVR+DPVY*PVY+DPRY*PRY+DVRY*VRY+DPPP*PPP+DVVV*VVV 

+DRRR*RRR+DYYY*YYY 



PPV 

PVV 

PPR 

PRR 

PPY 

PYY 

WR 

VRR 

VVY 

VYY 

RRY 

RYY 

PVR 

PVY 

PRY 

VRY 

PPP 

VVV 

RRR 

YYY 



M12*M12 
Mll*M22 
M12*M12 
M11 *m 32 
M12*M12 
Mll*M42 
M22*M22 
M21*M32 
M22*M22 
M21*M42 
M32*M32 
M31*M42 
M12*M22 
M12*M22 
M12*M42 
M22*M42 
3.0*M11 
3.0 *m21 
3.0*M31 
3 . 0*M41 



*M21 

*M22 

*M31 

*M32 

*M41 

*M42 

*M31 

*M32 

*M41 

*M42 

*M41 

*M42 

*M31 

*M41 

*M31 

*M31 

*M12* 

*m22* 

*M32* 

*M42* 



M32 

M42 

M32 

M32 



M12 

M22 

M32 

M42 



★ Mil* 

*1412* 

*1411* 

*M12* 

*Mll* 

*M12* 

*M21* 

*M22* 

*M21* 

*M22* 

*M31* 

*M32* 

* (Mil 
*(M11 
*(M11 
*(M21 



m12*M22 

M21*M22 

M12*M32 

M31*M32 

M12*M42 

M41*M42 

M22*M32 

M31*M32 

M22*M42 

M41*M42 

M32*M42 

M41*M42 

*M22+M12*M21 ) 

*M22+M12*M21 ) 

*M42+M12*M41 ) 

*M42+M22*M41 ) 



L23=DPPV*PPV+DPVV*PVV+DPPR*PPR+DPRR*PRR+DPPY*PPY+DPYY*PYY 

+DVVR*VVR+DVRR*VRR+DVVY*VVY+DVYY*VYY+DRRY*RRY+DRYY*RYY 

+DPVR*PVR+DPVY*PVY+DPRY*PRY+DVRY*VRY+DPPP*PPP+DVVV*VVV 

+DRRR*RRR+DYYY*YYY 



L24= DPPV*M12*M12*M22 + DPVV*M12*M22*M22 
+ DPRR*M12*M32*M32 + DPPY*M12*M12*M42 
+ DVVR*M22*M22*M32 + DVRR*M22 *M32 *M3 2 
+ DVYY*M22*M42*M42 + DRRY*M32*M32*M42 
+ DPVR*M12*M22*m 32 + DPVY*Ml 2 *M2 2 *M4 2 
+ DVRY*M22*M42*M32 + DPPP*Ml2*Ml2*Ml2 
+ DRRR*M32*M32*M32 + DYYY*M42*M42*M42 



+ DPPR*M12*M12*M32 
+ DPYY*M12*M42*M42 
+ DVVY*M22*M22*M42 
+ DRYY*M32*M42*M42 
+ DPRY*M12*M42*M32 
+ DVVV*M22*M22*M22 



L31=L21 

L32«L22 

L33=L23 

L34=L24 

L21*L21*BB1*U*U 

L22=L22*BBl*U*U 

L23=L23*BB1*U*U 

L24«L24*BB1*U*U 

L31=L31*BB2*U*U 

L32*L32*BB2*U*U 
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L34«L34*BB2*U*U 

L41— 0 . (M21+U*Mll/3 . 0 ) 

L42=-M11* (M12*M21+0 . 5*Mll*M22+0 . 5*U*M12*M12) 
L43=-M12* (Mll*M22+0 . 5*Ml2*M21+0 . 5*U*M11*M12 ) 
L44«-0. 5*M12*M12* (M22+U*M12/3 . 0 ) 

Invert Transformation Matrix 

DO 2 1=1,4 
DO 3 J=l,4 

TINV( I , J)=0.0 
TSAVE(I,J)=T(I,J) 

3 CONTINUE 
2 CONTINUE 

CALL DLUD (4,4, TSAVE , 4 , TLUD , I VLUD ) 

DO 4 1=1,4 

IF ( IVLUD( I ) . EQ. 0 ) STOP 3003 

4 CONTINUE 

CALL DILU( 4 , 4 ,TLUD, IVLUD, SVLUD) 

DO 8 1=1,4 
DO 9 J=l,4 

TINV( I , J )=TLUD( I , J ) 

9 CONTINUE 

8 CONTINUE 

Check Inv(T)*A*T 

IMULT=2 

IF (IMULT.EQ.l) CALL MULT ( TINV , ASAVE , T ) 

IF (IMULT.EQ.O) STOP 

Definition of Nij 

N11=TINV( 1 , 1 ) 

N21=TINV(2,1) 

N31=TINV(3,1) 

N41=TINV( 4,1) 

N12=TINV( 1,2) 

N22=TINV( 2,2) 

N32=TINV( 3,2) 

N42=TINV( 4,2) 

N13=TINV( 1,3) 

N23=TINV( 2,3) 

N33=TINV( 3 , 3 ) 

N43=TINV( 4,3) 

N14=TINV( 1,4) 

N24=TINV(2,4) 

N34=TINV( 3,4) 

N44=TINV( 4,4) 

R11=N12*L21+N13*L31+N14*L41 

R12=N12*L22+N13*L32+N14*L42 

R13=N12*L23+N13*L33+N14*L43 

R14=N12*L24+N13*L34+N14*L44 

R21=N22*L21+N23*L31+N24*L41 

R22=N22*L22+N23*L32+N24*L42 

R23=N22*L23+N23*L33+N24*L43 

R24=N22*L24+N23*L34+N24*L44 
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tivaiuate Aipna' ana umega' 



DINC-0 . 001 
XDR -XD+DINC 
XDL *XD-DINC 
XD =XDR 
PSI «0.D0 

Y =0.D0 

V =0.D0 
DPSI-KlP 
DV -K2P 
DR =K3 

DY -K1*XD/(XD**2+Y**2 ) 

A(1,1)-0.0D0 

A(1,2)-0.0D0 

A(1,3)=1.0D0 

A(1,4)*0.0D0 

A(2,l)= BB1*U*U*DPSI 

A( 2 , 2 )=AA11*U+BB1*U*U*DV 
A( 2, 3 )=AA12*U+BB1*U*U*DR 
A( 2, 4 )= BBl*U*U*DY 

A(3,l)= BB2*U*U*DPSI 

A( 3 , 2 )=AA21*U+BB2*U*U*DV 
A( 3 , 3 )=AA22*U+BB2*U*U*DR 
A(3,4)* BB2*U*U*DY 

A(4,l )=U*DCOS(PSI )-V*DSIN(PSI ) 

A(4,2)= DCOS(PSI) 

A(4,3)=0.0D0 
A(4,4 )=0.0D0 

CALL RG( 4 , 4 ,A,WR,WI , 0 , ZZZ , I Vl , FVl , lERR) 
CALL DSTABL(DEOS,WR,WI,FREQ) 

ALPHR=DEOS 

OMEGR*FREQ 

XD=XDL 
PSI =0.D0 

Y =0.D0 

V =0.D0 
DPSI=KlP 
DV =K2P 
DR =K3 

DY =K1*XD/(XD**2+Y**2 ) 

A(l,l )=0.0D0 
A(1,2)=O.ODO 
A(1,3)*1.0D0 
A(1,4)=0.0D0 

A(2,l)= BBl*U*U*DPSI 

A(2,2)=AA11*U+BB1*U*U*DV 
A(2, 3)=AA12*U+BB1*U*U*DR 
A(2,4)= BB1*U*U*DY 

A(3,l)= BB2*U*U*DPSI 

3. 2 ) =AA21*U+BB2*U*U*DV 

3. 3 ) »AA22*U+BB2*U*U*DR 

^(3,4)= BB2*U*U*DY 

4 , 1 )=U*DCOS( PSI )-V*DSIN( PSI ) 

\(4,2)= DCOS(PSI) 

\(4,3)=0.0D0 

\(4,4)=0.0D0 
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CALL DSTABL{ DECS, WR,WI, FREQ) 

ALPHL-DEOS 

OMEGL-FREQ 



DALPHA* ( ALPHR-ALPHL ) /( XDR-XDL ) 

^ DOMEGA= ( OMEGR-OMEGL ) / ( XDR-XDL ) 

Evaluation of Hopf Bifurcation Coefficients 

COEFl=3 . 0*Rll+Rl3+R22+3 . 0*R24 
COEF2=3 . 0*R21+R23-Rl2-3 . 0*Rl4 
AMU2 *-COEFl/(8.0*DALPHA) 

BETA2=0.25*COEF1 

TAU2 =-( COEF2-DOMEGA*COEFl/DALPHA)/( 8 . 0*OMEGA0 ) 

PERD =2 . 0*3 . 1415927/OMEGAO 

PER -PERD*U/L 

X=XD/L 

K4»K1/XD 

I WRITE (11,*) TC,COEFl 

t WRITE (12,*) TC,PER 

I WRITE (13,*) TC,DALPHA 

TC=TC+0 . 015 
40 CONTINUE 
C TL»TL+0.02 

-C 50 CONTINUE 



" 1001 FORMAT (' ENTER TIME CONSTANT') 

1005 FORMAT (' ENTER TIME LAG') 

1006 FORMAT (' ENTER DO') 

- 2002 FORMAT (5E15.5) 

: END 

SUBROUTINE DSOMEG ( UK , WR , WI , OMEGA , CHECK ) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(4),WI(4) 

CHECK=-1.0D+25 
DO 1 1=1,4 

IF (WR( I ) .LT.CHECK) GO TO 1 
CHECK=WR( I ) 

I J=I 



1 CONTINUE 

OMEGA=DABS(WI(IJ) ) 
IF (WI( IJ) .GT.O.DO) 
IF (WI ( I J ) . LT. 0 . DO ) 
RETURN 
END 



UK- 

UK 



= IJ 
= IJ-1 






SUBROUTINE DSTABL ( DEOS , WR , WI , OMEGA) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(4),WI(4) 

DEOS=-l . OD+20 
DO 1 1=1,4 

IF (WR( I ) .LT.DEOS) GO TO 1 
DEOS=WR( I ) 

I J=I 

CONTINUE 

OMEGA=WI(IJ) 

OMEGA=DABS ( OMEGA) 
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UBROUTINE NORMAL(T) 

MPLICIT DOUBLE PRECISION (A-H,0-Z) 

IMENSION T(4,4) ,TN0R(4,4) 

FAC=T(1,1)**2+T(1,2)**2 

F ( DABS(CFAC) .LE. ( 1 .D-10 ) ) STOP 4001 

NOR(1,1)=1.DO 

N0R(2,1)=(T(1,1)*T(2,1)+T(2,2)*T(1,2) )/CFAC 
’N0R(3,1)-(T(1,1)*T(3,1)+T( 3,2)*T(1,2) )/CFAC 
'NOR (4,1)=(T(1,1)*T(4,1)+T(4,2)*T(1,2))/CFAC 
NOR(1,2)-O.DO 

•N0R(2,2)*(T(2,2)*T(1,1)-T(2,1)*T(1,2) )/CFAC 
*NOR( 3,2)=(T(3,2)*T(1,1)-T(3,1)*T(1,2) )/CFAC 
'NOR (4,2)-(T(4,2)*T(l,l)-T(4,l)*T(l,2))/CFAC 
‘0 1 1 = 1,4 
DO 2 J=l,2 

T(I,J)=TNOR(I,J) 

CONTINUE 

ONTINUE 

.ETURN 

|;nd 

UBROUTINE MULT( TINV,A,T) 

MPLICIT DOUBLE PRECISION (A-H,0-Z) 

•IMENSION TINV( 4 ,4),A(4,4),T(4,4),Al(4,4),A2(4,4) 
•0 1 1 = 1,4 
DO 2 J=l,4 
A1(I, J)=0.D0 
A2( I, J)=0.D0 
CONTINUE 
:ONTINUE 
>0 3 1 = 1,4 
DO 4 J=l,4 
DO 5 K=l,4 

Al( I, J)=A( I,K)*T(K, J)+A1(I, J) 

CONTINUE 
CONTINUE 
:ONTINUE 
)0 6 1 = 1,4 
DO 7 J=l,4 
DO 8 K=l,4 

A2( I , J)=TINV( I ,K)*A1( K, J)+A2( I , J) 

CONTINUE 
CONTINUE 
:ONTINUE 
)0 11 1 = 1,4 

WRITE (*,101) (A( I , J) , J=1 , 4 ) 

:ONTINUE 
)0 12 1 = 1,4 

WRITE (*,101) (T(I, J) , J=l,4) 

:ONTINUE 
)0 10 1 = 1,4 

WRITE (*,101) (A2(I, J) ,J=1,4) 

:ONTINUE 

WRITE (*,101) A2(l,l) 

RETURN 

FORMAT (4E15.5) 

:nd 
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APPENDIX B 



PROGRAM SIMU.FOR (VEHICLE PATH SIMULATION) 
REAL Kl , K2 , K3 , L , NR , NV , NDRS , NDRB ,1Z, MASS , 

& NRDOT , NVDOT , Kl P , K2 P 

OPEN (20,FILE='SIMl.DAT' , STATUS* ' NEW ' ) 

OPEN ( 30 , FILE*' SIM2 . DAT' , STATUS* 'NEW' ) 



WEIGHT*435 . 0 
IZ *45.0 
L =7.3 

RHO =1.94 
G *32.2 

XG =0.0104 
MASS *WEIGHT/G 
U =5.0 
RATIO =0.0 
DO =0.4 
XG =XG/L 

MASS =MASS/( 0 . 5*RHO*L**3 ) 

IZ =IZ/( 0 . 5*RHO*L**5 ) 

YRDOT =-0.00178 
YVDOT =-0.03430 
YR =+0.01187 
YV =-0.03896 
YDRS =+0.02345 
YDRB =+0.02345 
NRDOT =-0.00047 
NVDOT *-0.00178 
NR =-0.01022 
NV =-0.00769 
NDRS =-0.337*0.02345 
NDRB *+0 .283*0.02345 

DH =( IZ-NRDOT) * (MASS- YVDOT )- 
& (MASS*XG-YRDOT) * (MASS*XG-NVDOT) 

AA11=( ( IZ-NRDOT) *YV-( MASS *XG- YRDOT ) *NV)/DH 
AA12=( ( IZ-NRDOT) *(-MASS+YR)- 
& (MASS*XG-YRDOT) * (-MASS*XG+NR) ) /DH 

AA21=( (MASS-YVDOT) *NV-(MASS*XG-NVDOT) *YV)/DH 
AA22= ( ( MASS-YVDOT ) * ( -MASS*XG+NR ) - 
& (MASS*XG-NVDOT) * (-MASS+YR) )/DH 

BBll*( ( IZ-NRDOT) *YDRS-(MASS*XG-YRDOT) *NDRS)/DH 
BB12*( ( IZ-NRDOT) * YDRB- (MASS *XG- YRDOT ) *NDRB)/DH 
BB21= ( ( MASS-YVDOT ) *NDRS- ( MASS *XG-NVDOT ) * YDRS ) /DH 
BB22=( (MASS-YVDOT) *NDRB-( MASS *XG-NVDOT) * YDRB) /DH 

BBl=BBll+RATIO*BBl2 

BB2=BB21+RATIO*BB22 

WRITE (*,*) 'ENTER TC,TL,D' 
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)LE=1 . 0/TC 

n = 3.0*POLE 

)2-3.0*POLE**2 

)3=P0LE**3 

L=BB1 

L-BB2 

L*-ADl-(AAll+AA22) 
2*(BB1*AA22-BB2*AA12 ) 
2=(BB2*AAll-BBl*AA21 ) 
L=AD3/(BB2*AA11-BB1*AA21 ) 
2»AD2-(AA11*AA22-AA12*AA21 )+BB2*K1 
2-(C1*B2-C2*B1)/(A1*B2-A2 *b1 ) 
3«(C2*A1-C1*A2)/(A1*B2-A2*B1) 

IME=100 . 0 

r=o.i 

-TIME/DT 

Y=TL/DT 

SI= 0.0 

0.0 

= 0.0 
= 0.5 
= 0.0 
COUNT=0 
CON=Y 

0 10, 1=1, N 

T= T+DT 

PSIDOT=R 

VDOT =AA11*V + AA12*R + BBl*DEL 
RDOT =AA21*V + AA22*R + BB2*DEL 
XDOT =COS ( PSI )-V*SIN( PSI ) 

YDOT =SIN(PSI )+V*COS( PSI ) 

PSI=PSI+( PSIDOT*DT) 

V=V+( VDOT*DT) 

R=R+(RDOT*DT) 

X=X+(XDOT*DT) 

Y=Y+( YDOT*DT) 

ICOUNT=ICOUNT+l 

IF ( ICOUNT.NE. lY) GO TO 11 

YCON=Y 

ICOUNT=0 

PSIC=-ATAN( YCON/D) 

DEL=K1* ( PSI-PSIC)+K2*V+K3*R 
DEL=D0*TANH( DEL/DO ) 

WRITE(20,*) T,Y 
WRITE(30,*) T,Y 
CONTINUE 
END 
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APPENDIX C 



PROGRAM THESNEWT.FOR (TWO POINT TIME DELAY ANALYSIS) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION Kl , K2 , K3 , L , NR , NV , NDRS , NDRB , I Z , MASS , 
& NRDOT,NVDOT 

DIMENSION A( 4 , 4 ) , FVl ( 4 ) , IVl ( 4 ) , ZZZ( 4 , 4 ) , WR ( 4 ) , WI( 4 ) 
DIMENSION DIST(10,10) 

OPEN ( 11 , FILE='TCN8 .DAT' , STATUS= ' NEW ' ) 

Vehicle Parameters: 

U= 5.0 

RATIO= 0.0 
WEIGHT=435 . 0 
IZ =45.0 
L =7.3 

RHO =1.94 
G =32.2 
XG =0.0104 
MASS =WEIGHT/G 

YRDOT =-0 . 00178*0 . 5*RHO*L**4 
YVDOT =-0 . 03430*0 . 5*RHO*L**3 
YR =+0.01187*0. 5*RHO*L**3 
YV =-0 . 03896*0 . 5*RHO*L**2 
YDRS =+0.02345*0. 5*RHO*L**2 
YDRB =+0 . 02345*0 . 5*RHO*L**2 
NRDOT =-0 . 00047*0 . 5*RHO*L**5 
NVDOT =-0 . 00178*0 . 5*RHO*L**4 
NR =-0 . 01022*0 . 5*RHO*L**4 
NV =-0 . 00769*0 . 5*RHO*L** 3 
NDRS =-0.3 37*0. 0 2 3 4 5 * 0 . 5 *RHO* L* * 3 
NDRB =+0.283*0. 02345*0 . 5*RHO*L** 3 

DH = ( I Z-NRDOT ) * ( MASS-YVDOT ) - 
& (MASS*XG-YRDOT) * (MASS*XG-NVDOT) 

AA11=( ( IZ-NRDOT) *YV-(MASS*XG-YRDOT) *NV)/DH 
AA12=( ( IZ-NRDOT) * (-MASS+YR )- 
& (MASS*XG-YRDOT) * (-MASS*XG+NR) )/DH 

AA21=( (MASS-YVDOT) *NV-(MASS*XG-NVDOT) *YV)/DH 
AA22=( (MASS-YVDOT) *(-MASS*XG+NR)- 
& (MASS*XG-NVDOT)*(-MASS+YR) )/DH 

BB11=( ( IZ-NRDOT) * YDRS- (MASS *XG- YRDOT ) *NDRS )/DH 
BB12=( ( IZ-NRDOT) * YDRB- (MASS*XG-YRDOT) *NDRB)/DH 
BB21=( (MASS-YVDOT) *NDRS-(MASS*XG-NVDOT) *YDRS )/DH 
BB22= ( ( MASS-YVDOT ) *NDRB- ( MASS *XG-NVDOT ) * YDRB ) /DH 

BB1=BB11+RATI0*BB12 

BB2=BB21+RATIO*BB22 

WRITE (*,1005) 



68 



READ ( * , * ) 



TL 



RITE (*,1006) 

EAD ( * , * ) TC 

PS-1 .D-10 
TMAX-100000 

L-0.01 
0 4 1-1,100 
LD - TL * L/U 

DO 1 J-1,100 
TCD-TC*L/U 
POLEC-1 . 0/TCD 
AD1=3.0*POLEC 
AD2-3 . 0*POLEC**2 
AD3-POLEC**3 
Al=BBl*U*U 
B1-BB2*U*U 

Cl--ADl-( AA11+AA22 ) *U 
A2-(BB1*AA22-BB2*AA12 ) *U**3 
B2=(BB2*AA11-BB1*AA21 ) *U**3 
Kl-AD3/( ( BB2*AA11-BB1*AA21 ) *U**3 ) 

C2=AD2-( AA11*AA22-AA12*AA21 ) *U* *2+BB2*U*U*Kl 
K2=(C1*B2~C2*B1 )/(Al*B2-A2*Bl ) 

K3=(C2*A1-C1*A2 )/(Al*B2-A2*Bl ) 

A1- (AA11*BB2-AA21*BB1 ) *U*U*U*K1 

A2- (AA11*AA22-AA12*AA21 ) *U*U+ ( AAl 1 *BB2-AA21 *BBl ) *U*U*U*K3+ 
(AA22*BB1-AA12*BB2 ) *U*U*U*K2-BB2*U*U*K1 
A3— (AA11+AA22 ) *U- ( BBl *K2+BB2*K3 ) *U*U 
BO- (AA11*BB2-AA21*BB1 )*U*U*U*U*Kl 
B1=-(AA12*BB2-AA22*BB1 ) *U*U*U* K1-BB2 *U*U*U* Kl 
B2 — BB1*U*U*K1 

AQ-B2*A3-B1 

BQ=B1*A2-B2*A1-B0*A3 

CQ=B0*A1 

DET-BQ**2-4 .D0*AQ*CQ 
IF (DET.LT.O.DO) GO TO 4 
ZT1=(-BQ+DSQRT(DET) )/( 2 . 0D0*AQ) 

ZT2- ( -BQ-DSQRT ( DET ) )/( 2 . ODO *AQ ) 

IF (ZTl.GE.O.DO) 0M=DSQRT(ZT1) 

IF (ZT2.GE.0.D0) OM-DSQRT ( ZT2 ) 

ALPHAl-AQ 

ALPHA2-BQ 

ALPHA3-CQ 

BETA1--B2 

BETA2-B0+B2*A2-B1*A3 

BETA3-Bl*Al-B0*A2 

OMOLD-OM 

BETA-BETAl*OMOLD**5+BETA2*OMOLD**3+BETA3*OMOLD 
ALPHA-ALPHA1*0M0LD**4+ALPHA2*0M0LD**2+ALPHA3 
BEP-5. D0*BETAl*OMOLD**4+3 .DO*BETA2*OMOLD**2+BETA3 
ALP-4 . DO*ALPHAl*OMOLD**3+2.DO*ALPHA2*OMOLD 
F1-BETA*DSIN( 2.DO*OMOLD*TLD)-2.DO*BETA*DSIN(OMOLD*TLD) 
F2-ALPHA*DC0S ( 2 . DO*OMOLD*TLD ) -2 . DO *ALPHA*DCOS ( OMOLD*TLD ) 
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FP1-BEP*DSIN( 2 ,DO*OMOLD*TLD)+BETA*2 . DO*TLD*DCOS ( 2 . DO*OMOLD*TLD ) 

FP2»-2 .DO*BEP*DSIN(OMOLD*TLD)-2 .DO*BETA*TLD*DCOS(OMOLD*TLD) 

FP3«ALP*DC0S( 2 ,DO*OMOLD*TLD)-ALPHA*2 . DO*TLD*DSIN( 2 . DO*OMOLD*TLD ) 

FP4=-2 ,DO*ALP*DCOS(OMOLD*TLD)+2.DO*TLD*ALPHA*DSIN(OMOLD*TLD) 

FPRIME» FP1+FP2+FP3+FP4 

IF (FPRIME.EQ.O.DO) STOP 1112 

IF (F.EQ.O.DO) GO TO 2 

OMNEW-OMOLD-F/FPRIME 

OMDI F=DABS ( OMNEW-OMOLD ) 

IF (OMDIF.LE.EPS) GO TO 2 

OMOLD=OMNEW 

IT=IT+1 

IF ( IT.GT. ITMAX) STOP 1111 
GO TO 3 

2 OM=OMNEW 

XDNUMl=DSQRT( ( B0-B2 *OM* * 2 ) * * 2+ ( Bl *OM ) * * 2 ) 

XDNUM2=DSQRT ( ( 2 . DO *DCOS ( OM*TLD ) -DCOS ( 2 . DO *OM*TLD ) ) * * 2 + 

& (DSIN( 2 ,DO*OM*TLD)-2 .DO*DSIN(OM*TLD) ) **2 ) 

' XDNUM»XDNUMl *XDNUM2 

XDDEN=OM*DSQRT( ( A1-A3 *OM* * 2 ) * * 2+ ( A2 *OM-OM* * 3 ) * * 2 ) 

XD=XDNUM/XDDEN 
X=XD/L 
DIST(I,J)=X 
WRITEdl,*) X,TL 
: TC=TC + 0.015 

: 1 CONTINUE 

TL=TL + 0.02 
4 CONTINUE 

: WRITE (11, 20) ( (DIST(I,J) ,J=1,10) ,1=1,10) 

20 FORMAT( 10 ( 2X, F4 . 2 ) ) 

1005 FORMAT [' ENTER TIME LAG') 

1006 FORMAT (' ENTER TIME CONSTANT') 

END 
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